TempOrder_Exp3_Analysis

Dataset: 64 participants collected in 2024/03 on SONA. 15 Excluded, totaling 49 used in the analysis.

# options(conflicts.policy = "strict")
# Basic
library(dplyr, mask.ok = c("filter"))
#library(MASS)
library(tidyverse)
# Analysis
library(lme4)
library(lmerTest)
library(performance)
library(jtools)
library(optimx)
library(lsr)
library(emmeans)
library(car)  
library(brglm2)
library(multcomp)
# Plotting
library(ggplot2)
library(lattice) 
library(sjPlot)
library(ggpubr)
library(plotrix)
library(DHARMa)
library(brms)

# Disable scientific notation
options(scipen=999) 

theme_set(theme_classic())

# GLMER Optimizers
cl1 <- glmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl = list(maxfun = 1e+09), check.conv.grad = .makeCC("warning",
tol = 0.002, relTol = NULL), check.conv.singular = .makeCC(action = "message", tol = 1e-09), check.conv.hess = .makeCC(action = "warning", tol = 0.000001))
cl2 <- glmerControl(optimizer = "Nelder_Mead", calc.derivs = FALSE, optCtrl = list(maxfun = 1e+09), check.conv.grad = .makeCC("warning", tol = 0.002, relTol = NULL), check.conv.singular = .makeCC(action = "message", tol = 1e-09), check.conv.hess = .makeCC(action = "warning", tol = 0.000001))
cl3 <- glmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "nlminb", maxiter = 1e+09), check.conv.grad = .makeCC("warning", tol = 0.002, relTol = NULL), check.conv.singular = .makeCC(action = "message", tol = 1e-09), check.conv.hess = .makeCC(action = "warning", tol = 0.000001))
cl4 <- glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE, optCtrl = list(maxfun = 1e+09), check.conv.grad = .makeCC("warning", tol = 0.002, relTol = NULL), check.conv.singular = .makeCC(action = "message",
tol = 1e-09), check.conv.hess = .makeCC(action = "warning",
tol = 0.000001))

# LMER Optimizers
cl_a <- lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+09))
cl_b <- lmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 1e+09))
cl_c <- lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb",
    maxiter = 1e+09))
cl_d <- lmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 1e+09))

1. Data Prep

1.1 Import data files

# set the path of the folder containing CSV files
folder_path <- "/Users/dcl/Library/CloudStorage/Box-Box/DCL_ARCHIVE/Documents/Events/exp165_TemporalOrder/Manuscript/Complete_Data/Data_Deidentification/Data_Exp3"

# get the list of CSV files in the folder
csv_files <- list.files(folder_path, pattern = "*.csv")

# initialize an empty list to store dataframes
df_list <- list()

# loop through the list of CSV files, reading each file into a dataframe and adding it to the list
for (csv_file in csv_files) {
  df <- read.csv(file.path(folder_path, csv_file))
  df_list[[csv_file]] <- df
}

# combine all dataframes in the list into a single dataframe
combined_df <- do.call(rbind, df_list)

demographics_combined_df <- combined_df %>% filter(trial_index == 2) %>% 
  dplyr::select('run_id', 'response')

# Helper Function 1: Extract Age
extract_text_between_age_and_gender <- function(text) {
  pattern <- "(?<=\\bage\\b)(.*?)(?=\\bgender\\b)"
  matches <- regmatches(text, gregexpr(pattern, text, perl = TRUE))
  unlist(matches)
}
# Helper Function 2: Extract Gender
extract_text_between_gender_and_race <- function(text) {
  pattern <- "(?<=\\bgender\\b)(.*?)(?=\\brace\\b)"
  matches <- regmatches(text, gregexpr(pattern, text, perl = TRUE))
  unlist(matches)
}
# Helper Function 3: Eliminate Text 
eliminate_text <- function(text){
  gsub("\\D", "", text)
}
# Helper Function 4: Eliminate Symbols
eliminate_symbols <- function(text){
  gsub("[[:punct:]]", "", text)
}

# Helper Function 5: Count "Yes"
count_yes <- function(input_string) {
  matches <- gregexpr("Yes", input_string, ignore.case = TRUE)
  match_count <- sum(sapply(matches, function(x) if (x[1] == -1) 0 else length(x)))
  return(match_count)
}
  
# Helper Function 6: Extract Race
extract_text_between_race_and_hispLat <- function(text) {
  pattern <- "(?<=\\brace\\b)(.*?)(?=\\bhispanicLatino\\b)"
  matches <- regmatches(text, gregexpr(pattern, text, perl = TRUE))
  unlist(matches)
}

# Helper Function 7: Extract Education_years
extract_education_years <- function(text) {
  return(substr(text, 21, 22))
}
  
# Helper Function 8: Extract HispanicLatino
extract_hispLat <- function(text) {
  string_length <- nchar(text)
  return(substr(text, string_length - 5, string_length - 1))
}

# Get Age Information
demographics_combined_df$age_raw <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_text_between_age_and_gender)
demographics_combined_df$age <- apply(demographics_combined_df[, "age_raw", drop = FALSE], 1, eliminate_text)
demographics_combined_df$age <- as.numeric(demographics_combined_df$age)

# Get Gender Information
demographics_combined_df$gender_raw <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_text_between_gender_and_race)
demographics_combined_df$gender <- apply(demographics_combined_df[, "gender_raw", drop = FALSE], 1, eliminate_symbols)

# Get Race Information
demographics_combined_df$race_raw <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_text_between_race_and_hispLat)
demographics_combined_df$race <- apply(demographics_combined_df[, "race_raw", drop = FALSE], 1, eliminate_symbols)

# Get Education_years Information
demographics_combined_df$edu_years <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_education_years)

# Get Hispanic/Latino Information
demographics_combined_df$hispLat_raw <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_hispLat)
demographics_combined_df$hispLat <- apply(demographics_combined_df[, "hispLat_raw", drop = FALSE], 1, eliminate_symbols)

# Select only useful columns
demographics_df <- demographics_combined_df %>% 
  dplyr::select('run_id', 'age', 'gender', 'race', 'edu_years', 'hispLat')

# Get Gender Distribution
gender_df <- as.data.frame(table(demographics_df$gender))
gender_df <- gender_df %>% rename(Gender = Var1, 
                                  Number = Freq)
gender_df$percentage <- gender_df$Number / nrow(demographics_df)
gender_df

# Get Race Distribution
race_df <- as.data.frame(table(demographics_df$race))
race_df <- race_df %>% rename(Race = Var1, 
                              Number = Freq)
race_df$percentage <- race_df$Number / nrow(demographics_df)
race_df

# Get Hispanic/Latino Distribution
hispLat_df <- as.data.frame(table(demographics_df$hispLat))
hispLat_df <- hispLat_df %>% rename(HispanicLatino = Var1, 
                                    Number = Freq)
hispLat_df$percentage <- hispLat_df$Number / nrow(demographics_df)
hispLat_df

edu_years_df <- as.data.frame(table(demographics_df$edu_years))
edu_years_df <- edu_years_df %>% rename(edu_years = Var1, 
                                    Number = Freq)
edu_years_df$percentage <- edu_years_df$Number / nrow(demographics_df)
edu_years_df

mean(as.numeric(demographics_df$edu_years))

## [1] 13.34375

sd(as.numeric(demographics_df$edu_years))

## [1] 1.311957

mean(demographics_combined_df$age)

## [1] 19.71875

sd(demographics_combined_df$age)

## [1] 1.174717

nrow(demographics_combined_df)

## [1] 64

min(demographics_combined_df$age)

## [1] 18

max(demographics_combined_df$age)

## [1] 23

demographics_combined_df$gender_raw <- apply(demographics_combined_df[, "response", drop = FALSE], 1, extract_text_between_gender_and_race)
demographics_combined_df$gender <- apply(demographics_combined_df[, "gender_raw", drop = FALSE], 1, eliminate_symbols)

table(demographics_combined_df$gender)

## 
## Female   Male 
##     46     18

complete_data <- combined_df %>% dplyr::select("run_id", "response", "rt", "time_elapsed", "stimuli_order", "trial_index", "task", "stimuli_type", "stimuli_name", "RQ_question_index", "RQ_correct_answer", "pair_num", "pair_location", "pair_first", "pair_second", "pair_type", "normal0_shuffle1", "correct_answer", "SourceMem_question_index", "pair_within_index", "coarse_count", "SourceMem_correct_answer", "Source_Sentence_index", "SourceMem_correct", "prompt", "choices", "answer_clean", "answer_source" )

2. Data Exclusion

2.1 Removing individual participants

Criteria 1: Exclude - Subjects who reported technical problems (Excluded = 0)

#complete_data %>% 
#  filter(trial_index == 645) %>% 
#  dplyr::select(run_id, response)

# Store problematic IDs for Criteria 1 (reported technical problems)
exclusion_c1_list <- c()

Criteria 2: Exclude - Subjects whose reaction time > 40000 ms for more than 5 individual trials or < 300 ms for more than 5 individual trials that require meaningful responses

# Frequency Histogram of RT distribution prior to data exclusion
complete_data$rt <- as.numeric(complete_data$rt)
complete_data %>%
  filter(task == 'encoding_sentence' | task == 'TOM' | task == 'TOM_Conf' | task == 'TOM_Coarse' | task == 'TOM_Coarse_Conf'| task == 'SourceMem') %>%
  ggplot(aes(x=rt)) +
    geom_histogram(binwidth = 500) +
    coord_cartesian(xlim = c(0, 70000)) +
    labs(title = "Frequency Histogram of RT distribution prior to data exclusion", x = "Reaction Time", y = "Count")

# Filter out problematic IDs (problematic: rt > 40000 ms for > 5 trials, or < 300 ms for > 5 trials)
trial_too_long <- complete_data %>% filter(task == 'encoding_sentence' | task == 'TOM' | task == 'TOM_Conf' | task == 'TOM_Coarse' | task == 'TOM_Coarse_Conf'| task == 'SourceMem') %>% filter(rt > 40000) %>% group_by(run_id) %>% summarize(num_trial_too_long = sum(rt > 40000)) %>% filter(num_trial_too_long > 5)
trial_too_long 

trial_too_short <- complete_data %>% filter(task == 'encoding_sentence' | task == 'TOM' | task == 'TOM_Conf' | task == 'TOM_Coarse' | task == 'TOM_Coarse_Conf'| task == 'SourceMem')%>% filter(rt < 300) %>% group_by(run_id) %>% summarize(num_trial_too_short = sum(rt < 300)) %>% filter(num_trial_too_short > 5)
trial_too_short 

# Store problematic IDs for Criteria 2.1 (rt > 40000 ms for > 5 trials)
exclusion_c2.1_list <- trial_too_long$run_id 

# Store problematic IDs for Criteria 2.2 (rt < 300 ms for > 5 trials)
exclusion_c2.2_list <- trial_too_short$run_id  

Criteria 3: Exclude subjects whose mean accuracy for RQ Questions < 60%

RQ_data <- complete_data %>% 
  filter(task == 'RQ') 

# use grepl() to check if "response" contains "correct_answer"
# use rowwise() to apply to each row
RQ_data <- RQ_data %>%
  rowwise() %>%
  mutate(RQ_response = grepl(RQ_correct_answer, response, ignore.case = TRUE)) %>% 
  mutate(RQ_correct = ifelse(RQ_response == "TRUE", 1, 0))

RQ_data$RQ_question_index_sep <- paste(RQ_data$stimuli_name, RQ_data$RQ_question_index, sep = "_")

# Frequency Histogram of Math Accuracy distribution prior to data exclusion - By subject
RQ_data_mean_bySub <- RQ_data %>% 
  group_by(run_id) %>% 
  summarize(meanRQ = mean(RQ_correct)) %>% 
  arrange(desc(meanRQ))

# RQ_data_mean_bySub # Show subject level grouping

RQ_data_mean_bySub %>% 
  ggplot(aes(x = meanRQ)) +
  geom_histogram() + 
  labs(title = "Frequency Histogram of Reading Comprehension Question Accuracy Distribution (by subject)", x = "Mean RQ Accuracy", y = "count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Filter out subjects that have < 60% mean accuracy for RQ questions
RQ_too_low <- RQ_data_mean_bySub %>% filter(meanRQ < 0.6) 

# Store problematic IDs for Criteria 3 (Mean RQ Accuracy < 60%)
exclusion_c3_list <- RQ_too_low$run_id

Criteria 4.1: Exclude - Subjects whose Fine-level Temporal Order Memory Mean Accuracy is < 50%

## Fine_TOM_raw_data: only Fine_TOM data
Fine_TOM_raw_data <- complete_data %>% 
  filter(task == 'TOM')

# Scoring: compare response with correct answer to get Fine_TOM (0/1)
Fine_TOM_raw_data$TOM_response <- as.numeric(Fine_TOM_raw_data$response)
Fine_TOM_raw_data$TOM_result <- as.numeric(Fine_TOM_raw_data$response == Fine_TOM_raw_data$correct_answer)

# Calculate TOM accuracy
Fine_TOM_raw_data_bySub <- Fine_TOM_raw_data %>% group_by(run_id) %>% summarize(meanTOM = mean(TOM_result))
Fine_TOM_raw_data_bySub <- ungroup(Fine_TOM_raw_data_bySub)
Fine_TOM_raw_data_bySub <- as.data.frame(Fine_TOM_raw_data_bySub)

# Frequency Histogram of mean TOM accuracy prior to data exclusion
Fine_TOM_raw_data_bySub %>% ggplot(aes(x=meanTOM)) +
    geom_histogram() + 
  labs(title = "Frequency Histogram of Mean Fine_TOM Accuracy prior to data exclusion", x = "Mean TOM Accuracy", y = "Count")

# Lower Limit: Exclude no one
trial_Fine_TOM_too_low <- Fine_TOM_raw_data_bySub %>% filter(meanTOM < 0.5) 

# Store problematic IDs for Criteria 4.1 (Mean TOM Accuracy < 0.5)
exclusion_c4.1_list <- trial_Fine_TOM_too_low$run_id 

#Fine_TOM_raw_data_bySub %>% arrange(desc(meanTOM))

Criteria 4.2: Exclude - Subjects whose Coarse-level Temporal Order Memory Mean Accuracy is < 50%

## Coarse_TOM_raw_data: only Coarse_TOM data
Coarse_TOM_raw_data <- complete_data %>% 
  filter(task == 'TOM_Coarse')

# Scoring: compare response with correct answer to get Fine_TOM (0/1)
Coarse_TOM_raw_data$TOM_response <- as.numeric(Coarse_TOM_raw_data$response)
Coarse_TOM_raw_data$TOM_result <- as.numeric(Coarse_TOM_raw_data$response == Coarse_TOM_raw_data$correct_answer)

# Calculate TOM accuracy
Coarse_TOM_raw_data_bySub <- Coarse_TOM_raw_data %>% group_by(run_id) %>% summarize(meanTOM = mean(TOM_result))
Coarse_TOM_raw_data_bySub <- ungroup(Coarse_TOM_raw_data_bySub)
Coarse_TOM_raw_data_bySub <- as.data.frame(Coarse_TOM_raw_data_bySub)

# Frequency Histogram of mean TOM accuracy prior to data exclusion
Coarse_TOM_raw_data_bySub %>% ggplot(aes(x=meanTOM)) +
    geom_histogram() + 
  labs(title = "Frequency Histogram of Mean Coarse_TOM Accuracy prior to data exclusion", x = "Mean TOM Accuracy", y = "Count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Lower Limit: Exclude no one
trial_Coarse_TOM_too_low <- Coarse_TOM_raw_data_bySub %>% filter(meanTOM < 0.5)

# Store problematic IDs for Criteria 4.2 (Mean Coarse_TOM Accuracy < 0.5)
exclusion_c4.2_list <- trial_Coarse_TOM_too_low$run_id 

#Coarse_TOM_raw_data_bySub %>% arrange(desc(meanTOM))

Criteria 4.3: Exclude - Subjects whose Source Memory Mean Accuracy is < 50%

## SourceMem_raw_data: only Coarse_TOM data
SourceMem_raw_data <- complete_data %>% 
  filter(task == 'SourceMem') %>% 
  rowwise() %>%
  # score source mem
  mutate(source_mem_TF =  grepl(SourceMem_correct_answer, response)) %>% 
  mutate(source_mem_result = ifelse(source_mem_TF == "TRUE", 1, 0))

# Calculate TOM accuracy
SourceMem_raw_data_bySub <- SourceMem_raw_data %>% group_by(run_id) %>% summarize(meanSourceMem = mean(source_mem_result))
SourceMem_raw_data_bySub <- ungroup(SourceMem_raw_data_bySub)
SourceMem_raw_data_bySub <- as.data.frame(SourceMem_raw_data_bySub)

# Frequency Histogram of mean TOM accuracy prior to data exclusion
SourceMem_raw_data_bySub %>% ggplot(aes(x = meanSourceMem)) +
    geom_histogram() + 
  labs(title = "Frequency Histogram of Mean Source Memory Accuracy prior to data exclusion", x = "Mean Source Memory Accuracy", y = "Count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Lower Limit
trial_SourceMem_too_low <- SourceMem_raw_data_bySub %>% filter(meanSourceMem < 0.5) 

# Store problematic IDs for Criteria 4.3 (Mean SourceMem Accuracy < 0.5)
exclusion_c4.3_list <- trial_SourceMem_too_low$run_id 

#SourceMem_raw_data_bySub %>% arrange(desc(meanSourceMem))

# Correlation between mean RQ and Fine_TOM accuracy (by Subject)
RQ_Fine_TOM_bySubject <- merge(RQ_data_mean_bySub, Fine_TOM_raw_data_bySub, by = "run_id")
RQ_Fine_TOM_bySubject %>% 
  ggplot(aes(x = meanRQ, y = meanTOM)) +
  geom_point() + 
  labs(title = "Correlation between mean RQ and Fine_TOM accuracy (by Subject)")

# Correlation between mean RQ and Coarse_TOM accuracy (by Subject)
RQ_Coarse_TOM_bySubject <- merge(RQ_data_mean_bySub, Coarse_TOM_raw_data_bySub, by = "run_id")
RQ_Coarse_TOM_bySubject %>% 
  ggplot(aes(x = meanRQ, y = meanTOM)) +
  geom_point() + 
  labs(title = "Correlation between mean RQ and Coarse_TOM accuracy (by Subject)")

# Correlation between mean RQ and Coarse_TOM accuracy (by Subject)
RQ_SourceMem_bySubject <- merge(RQ_data_mean_bySub, SourceMem_raw_data_bySub, by = "run_id")
RQ_SourceMem_bySubject %>% 
  ggplot(aes(x = meanRQ, y = meanSourceMem)) +
  geom_point() + 
  labs(title = "Correlation between mean RQ and Source Memory accuracy (by Subject)")

Criteria 5: Exclude - Subjects who had a delay time longer than 35 min

# Time elapsed when the 10th story of the encoding phase ended 
time_at_1st_test <- complete_data %>% subset(trial_index == 315) 
# Time elapsed when the 1st story of the encoding phase ended
time_at_1st_story <- complete_data %>% subset(trial_index == 35)

# Record each participant's delay time
time_diff <- as.data.frame(time_at_1st_test$time_elapsed - time_at_1st_story$time_elapsed)
time_diff <- as.data.frame(cbind(time_at_1st_test$run_id, time_diff$`time_at_1st_test$time_elapsed - time_at_1st_story$time_elapsed`))
# convert delay time to minutes
time_diff$diff <- time_diff$V2 / 1000 / 60


# display the distribution of delay time
time_diff %>% ggplot(aes(diff)) +
  geom_histogram() +
  labs(title = "Frequency Histogram of Delay Time prior to data exclusion", x = "Delay Time", y = "Count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

time_diff %>% arrange(desc(diff))

# Lower Limit: Exclude no one
trial_Delay_too_long <- time_diff %>% filter(diff > 35) 

# Store problematic IDs for Criteria 5 (Delay time > 35 min)
exclusion_c5_list <- trial_Delay_too_long$V1

Data Exclusion Decision:

# Decision: Merge ID lists from Criteria 1-4
exclusion_complete_list <- unique(c(exclusion_c1_list, exclusion_c2.1_list, exclusion_c2.2_list, exclusion_c3_list, exclusion_c4.1_list, exclusion_c4.2_list, exclusion_c4.3_list, exclusion_c5_list))
exclusion_complete_list # Print out IDs to be removed

##  [1]   6  12  14  39  44  54  80  87  88  96  98 119  50 111 120

# In "complete_data" dataframe, tag all the trials that are excluded based on subject exclusion criteria
complete_data$excluded_by_subject <- ifelse(complete_data$run_id %in% exclusion_complete_list, yes = "TRUE", no = "FALSE")

# Create "main_data" dataframe: Exclude problematic subjects, and filter to only include trials that are relevant in the analysis
main_data <- complete_data[!(complete_data$run_id %in% exclusion_complete_list), ] %>% filter(task == 'TOM' | task == 'TOM_Conf' | task == 'TOM_Conf' | task == 'TOM_Coarse' | task == 'TOM_Coarse_Conf' | task == 'SourceMem')
main_data$stimuli_type <- ifelse(test = main_data$stimuli_type == "HC", yes = "CS", no = "FS") # For stimuli_type, Change all "HC" to "CS" and all "LC" to "FS"

# Sample Size after data exclusion (n = 14)
nSample <- length(unique(main_data$run_id))
nSample # Print out sample size after removal

## [1] 49

# id_list after data exclusion
id_list <- unique(main_data$run_id)

across_trials <- main_data %>% filter(task == "TOM") %>% filter(pair_type == "across")

3. Data Preparation for Different Measures & Analyses for Source Memory and Recency Judgment

3.1 Source Memory Data

source_mem_data <- main_data %>% filter(task == "SourceMem")
# source_mem_result: Score if each source memory response is correct (correct = 1, incorrect = 0)
source_mem_data <- source_mem_data %>% 
  rowwise() %>%
  # score source mem
  mutate(source_mem_TF =  grepl(SourceMem_correct_answer, response)) %>% 
  mutate(source_mem_result = ifelse(source_mem_TF == "TRUE", 1, 0)) %>% 
  # create index
  mutate(order_mem_index = paste(stimuli_name, pair_num, sep = "_")) %>% 
  # person_trial_index: Make an index in the format of "runID_stimuliName_fineRecencyTrialIndex"
  # for merging data
  mutate(person_trial_index = paste(run_id, order_mem_index, sep = "_"))

Data exclusion: Exclude Source Memory trials that have reaction time < 300 ms, or > 3SD above the mean reaction time of all Source Memory trials

### Data exclusion
# Calculate sourceMem_rt_upperLim: 3SD above mean 
sourceMem_rt_upperLim <- mean(as.numeric(source_mem_data$rt)) + 3 * sd(as.numeric(source_mem_data$rt))
sourceMem_rt_upperLim  # Print out sourceMem_rt upper limit: 35277.93 ms

## [1] 35277.93

# Tag all the sourceMem trials that are excluded by rt < 300 or > 3SD above mean
source_mem_data$excluded_by_sourceMem_rt <- ifelse(as.numeric(source_mem_data$rt) > sourceMem_rt_upperLim | as.numeric(source_mem_data$rt) < 300, yes = TRUE, no = FALSE)
mean(source_mem_data$excluded_by_sourceMem_rt == TRUE) # Print out the proportion of sourceMem trials being excluded based on rt criterion

## [1] 0.005102041

3.1.1 Additional Analysis 1: Source Memory predicted by sentence location

** Additional Analysis 1**

Motivation: There are previous studies (Heusser et al., 2016) showing that source memory for boundary trials (at 1st position in a given event, after experiencing a perceptual boundary) are better than for non-boundary trials. Here, because each of our fine-level across-event pairs happens to contains sentence at the 3rd position in a coarse-level event and a sentence at the 1st position in the next adjacent pair, we have a chance to test if source memory accuracy is predicted by sentence location, controlling for narrative type (CS/FS).

Note that based on the nature of our stimuli, we believe that source memory in the current experiment is largely driven by participants’ semantic knowledge about given events, and therefore there should not be a significant main effect of sentence location during encoding. Since it is not the main focus of this experiment, we will not power our analysis based on this effect.

• Model: Mixed Effects Logistic Regression Model
• Outcome variable: Source Memory Accuracy (0/1)
• Predictors: stimuli_type (CS narrative vs. FS narrative) & within_coarse_loc (1st vs. 3rd)
• Random Effects: run_id (subject) + order_mem_index (different test pairs)

source_mem_data_analysis <- source_mem_data

source_mem_data_analysis$source_mem_result_binary <- as.factor(source_mem_data_analysis$source_mem_result)

# within_coarse_loc: whether the sentence is the 3rd sentence or the 1st sentence in its coarse-level event
source_mem_data_analysis <- source_mem_data_analysis %>% 
  mutate(within_coarse_loc = ifelse(pair_within_index == 1, 3, 1))

source_mem_data_analysis$within_coarse_loc <- as.factor(source_mem_data_analysis$within_coarse_loc)

# Descriptive Stats: 
# Mean Source Memory Accuracy by Sentence Location
source_mem_data_analysis %>% 
  group_by(within_coarse_loc) %>% 
  summarize(mean_sourceMem_accuracy = mean(source_mem_result),
            sd_source_Mem_accuracy = sd(source_mem_result))

# Mean Source Memory Accuracy by Sentence Location and Person
source_summ <- source_mem_data_analysis %>% 
  group_by(run_id, within_coarse_loc) %>% 
  summarize(mean_sourceMem_accuracy = mean(source_mem_result),
            sd_source_Mem_accuracy = sd(source_mem_result))

## `summarise()` has grouped output by 'run_id'. You can override using the
## `.groups` argument.

#source_summ
source_summ %>%
  ggplot(aes(x = within_coarse_loc, y = mean_sourceMem_accuracy, fill = within_coarse_loc)) +
    geom_violin() +
    geom_jitter(color="black", size=0.4, alpha=0.9, position = position_jitter(height = 0)) +
    geom_line(aes(group = run_id)) 

table(source_mem_data_analysis$stimuli_type, source_mem_data_analysis$source_mem_result_binary)

##     
##         0    1
##   CS  244 1716
##   FS  241 1719

mean(source_mem_data_analysis$source_mem_result)

## [1] 0.8762755

sd(source_mem_data_analysis$source_mem_result)

## [1] 0.329309

# 1. Select noInt vs. Int
# Model 1: With interaction between two predictors
source_mem_model_int <- glmer(source_mem_result_binary ~ stimuli_type * within_coarse_loc + (1 |run_id) + (1 + within_coarse_loc|order_mem_index), data = source_mem_data_analysis, family = binomial, control = cl1)

# Model 2: Without interaction between two predictors
source_mem_model_noInt <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 |run_id) + (1 + within_coarse_loc|order_mem_index), data = source_mem_data_analysis, family = binomial, control = cl1)

# Model selection favors Model 2
anova(source_mem_model_noInt, source_mem_model_int)

# Contrast Coding
source_mem_data_analysis_new <- source_mem_data_analysis
source_mem_data_analysis_new$within_coarse_loc <- as.factor(source_mem_data_analysis_new$within_coarse_loc)
source_mem_data_analysis_new$stimuli_type <- as.factor(source_mem_data_analysis_new$stimuli_type)

contrasts(source_mem_data_analysis_new$within_coarse_loc) <- contr.sum(2) 
contrasts(source_mem_data_analysis_new$stimuli_type) <- contr.sum(2)
contrasts(source_mem_data_analysis_new$within_coarse_loc)

##   [,1]
## 1    1
## 3   -1

contrasts(source_mem_data_analysis_new$stimuli_type)

##    [,1]
## CS    1
## FS   -1

# 2. Select random effect structure

# Model 1: Random Intercept only
source_mem_model_noInt_effCoding_a <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 |run_id) + (1|order_mem_index), data = source_mem_data_analysis_new, family = binomial, control = cl1)

# Model 2: Random Slope of within_coarse_loc on order_mem_index
source_mem_model_noInt_effCoding_b <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 |run_id) + (1 + within_coarse_loc|order_mem_index), data = source_mem_data_analysis_new, family = binomial, control = cl1)

# Model 3: Random Slope of within_coarse_loc on order_mem_index + Random Slope of stimuli_type on run_id
source_mem_model_noInt_effCoding_c <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 + stimuli_type |run_id) + (1 + within_coarse_loc|order_mem_index), data = source_mem_data_analysis_new, family = binomial, control = cl1)

# Model 4: Random Slope of within_coarse_loc on order_mem_index + Random Slope of within_coarse_loc on run_id
source_mem_model_noInt_effCoding_d <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 + within_coarse_loc|run_id) + (1 + within_coarse_loc|order_mem_index), data = source_mem_data_analysis_new, family = binomial, control = cl1)

# Model 5: All possible Random Slopes
source_mem_model_noInt_effCoding_e <- glmer(source_mem_result_binary ~ stimuli_type + within_coarse_loc + (1 + stimuli_type + within_coarse_loc |run_id) + (1 + stimuli_type + within_coarse_loc|order_mem_index), data = source_mem_data_analysis_new, family = binomial, control = cl1)

# Model selection favored model b
anova(source_mem_model_noInt_effCoding_a, source_mem_model_noInt_effCoding_b, source_mem_model_noInt_effCoding_c, source_mem_model_noInt_effCoding_d, source_mem_model_noInt_effCoding_e)

# Model Using Effect Coding (Ref Group: loc = 1, stimuli_type = CS)
source_mem_model <- source_mem_model_noInt_effCoding_b
summ(source_mem_model)

## MODEL INFO:
## Observations: 3920
## Dependent Variable: source_mem_result_binary
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 2101.47, BIC = 2145.39
## Pseudo-R² (fixed effects) = 0.01
## Pseudo-R² (total) = 0.56 
## 
## FIXED EFFECTS:
## -------------------------------------------------------
##                             Est.   S.E.   z val.      p
## ------------------------ ------- ------ -------- ------
## (Intercept)                 3.21   0.27    11.91   0.00
## stimuli_type1              -0.26   0.24    -1.10   0.27
## within_coarse_loc1          0.13   0.22     0.61   0.54
## -------------------------------------------------------
## 
## RANDOM EFFECTS:
## --------------------------------------------------
##       Group            Parameter        Std. Dev. 
## ----------------- -------------------- -----------
##      run_id           (Intercept)         0.80    
##  order_mem_index      (Intercept)         1.42    
##  order_mem_index   within_coarse_loc1     1.23    
## --------------------------------------------------
## 
## Grouping variables:
## -----------------------------------
##       Group        # groups   ICC  
## ----------------- ---------- ------
##      run_id           49      0.11 
##  order_mem_index      40      0.34 
## -----------------------------------

Anova(source_mem_model, type = "III")

# source_mem_data_s1: select all sourceMem trials of 1st sentence in each pair
source_mem_data_s1 <- source_mem_data %>%  
  filter(pair_within_index == 1)
# source_mem_data_s2: select all sourceMem trials of 2nd sentence in each pair
source_mem_data_s2 <- source_mem_data %>%  
  filter(pair_within_index == 2)

# change names of column "excluded_by_sourceMem_rt" to be s1/s2 specific
colnames(source_mem_data_s1)[colnames(source_mem_data_s1) == "excluded_by_sourceMem_rt"] <- "excluded_by_sourceMem1_rt"
colnames(source_mem_data_s2)[colnames(source_mem_data_s2) == "excluded_by_sourceMem_rt"] <- "excluded_by_sourceMem2_rt"

# source_mem_data_s2_selected: select only useful columns in source_mem_data_s2
source_mem_data_s2_selected <- source_mem_data_s2 %>% 
  dplyr::select(person_trial_index, source_mem_result, answer_source, coarse_count, excluded_by_sourceMem2_rt) 
# change column names in source_mem_data_s2_selected
colnames(source_mem_data_s2_selected) <- c("person_trial_index", "source_mem_result_s2", "answer_source_s2", "correct_source_s2", "excluded_by_sourceMem2_rt")
# MERGE source_mem_data_s2_selected with source_mem_data_s1
source_mem_data_merged <- merge(source_mem_data_s1, source_mem_data_s2_selected, by = "person_trial_index") 
# change column names for s1 columns in source_mem_data_merged
colnames(source_mem_data_merged)[colnames(source_mem_data_merged) == "source_mem_result"] <- "source_mem_result_s1"
colnames(source_mem_data_merged)[colnames(source_mem_data_merged) == "answer_source"] <- "answer_source_s1"
colnames(source_mem_data_merged)[colnames(source_mem_data_merged) == "coarse_count"] <- "correct_source_s1"

# source_mem_data_cleaned: select all useful columns for merging with other tests
source_mem_data_cleaned <- source_mem_data_merged %>% dplyr::select("person_trial_index", "source_mem_result_s1", "source_mem_result_s2", "answer_source_s1", "answer_source_s2", "correct_source_s1", "correct_source_s2", "excluded_by_sourceMem1_rt", "excluded_by_sourceMem2_rt")

3.2 Coarse-level Recency Judgment Data

Coarse_TOM_data <- main_data %>% filter(task == "TOM_Coarse")
Coarse_TOM_data <- Coarse_TOM_data %>% 
  # select only useful columns
  dplyr::select(-c("SourceMem_question_index", "pair_within_index", "coarse_count", "SourceMem_correct_answer", "Source_Sentence_index", "SourceMem_correct", "answer_source", "prompt", "choices", "answer_clean")) %>% 
  # create index
  mutate(order_mem_index = paste(stimuli_name, pair_num, sep = "_")) %>% 
  mutate(person_trial_index = paste(run_id, order_mem_index, sep = "_")) %>% 
  # score coarse-level TOM 
  mutate(Coarse_TOM_result = as.numeric(response == correct_answer))

Data exclusion: Exclude Coarse_TOM trials that have reaction time < 300 ms, or > 3SD above the mean reaction time of all Coarse_TOM trials

### Data exclusion
# Calculate Coarse_TOM_rt_upperLim: 3SD above mean 
Coarse_TOM_rt_upperLim <- mean(as.numeric(Coarse_TOM_data$rt)) + 3 * sd(as.numeric(Coarse_TOM_data$rt))
Coarse_TOM_rt_upperLim  # Print out sourceMem_rt upper limit: 17242.39 ms

## [1] 17242.39

# Tag all the Coarse_TOM trials that are excluded by rt < 300 or > 3SD above mean
Coarse_TOM_data$excluded_by_Coarse_TOM_rt <- ifelse(as.numeric(Coarse_TOM_data$rt) > Coarse_TOM_rt_upperLim | as.numeric(Coarse_TOM_data$rt) < 300, yes = TRUE, no = FALSE)
mean(Coarse_TOM_data$excluded_by_Coarse_TOM_rt == TRUE) # Print out the proportion of Coarse_TOM trials being excluded based on rt criterion

## [1] 0.01071429

# Coarse_TOM_data_cleaned: select all useful columns for merging with other tests
Coarse_TOM_data_cleaned <- Coarse_TOM_data %>% 
  dplyr::select("pair_first", "pair_second", "person_trial_index", "Coarse_TOM_result", "excluded_by_Coarse_TOM_rt")

## Coarse_TOM_Conf data
Coarse_TOM_Conf_data <- main_data %>% filter(task == "TOM_Coarse_Conf")
Coarse_TOM_Conf_data <- Coarse_TOM_Conf_data %>% 
  # select only useful columns
  dplyr::select(-c("SourceMem_question_index", "pair_within_index", "coarse_count", "SourceMem_correct_answer", "Source_Sentence_index", "SourceMem_correct", "answer_source", "prompt", "choices", "answer_clean")) %>% 
  # create index
  mutate(order_mem_index = paste(stimuli_name, pair_num, sep = "_")) %>% 
  mutate(person_trial_index = paste(run_id, order_mem_index, sep = "_")) 

Coarse_TOM_Conf_data_cleaned <- Coarse_TOM_Conf_data %>% dplyr::select("person_trial_index", "response")

# MERGE: Coarse_TOM + Coarse_TOM_Conf
Coarse_TOM_complete_data_cleaned <- merge(Coarse_TOM_data_cleaned, Coarse_TOM_Conf_data_cleaned, by = "person_trial_index") 

colnames(Coarse_TOM_complete_data_cleaned)[colnames(Coarse_TOM_complete_data_cleaned) == "response"] <- "Coarse_TOM_Conf"
colnames(Coarse_TOM_complete_data_cleaned)[colnames(Coarse_TOM_complete_data_cleaned) == "pair_first"] <- "Coarse_pair_first"
colnames(Coarse_TOM_complete_data_cleaned)[colnames(Coarse_TOM_complete_data_cleaned) == "pair_second"] <- "Coarse_pair_second"

3.3 Fine-level Recency Judgment Data

Fine_TOM_data <- main_data %>% filter(task == "TOM")
Fine_TOM_data <- Fine_TOM_data %>% 
  # select only useful columns
  dplyr::select(-c("SourceMem_question_index", "pair_within_index", "coarse_count", "SourceMem_correct_answer", "Source_Sentence_index", "SourceMem_correct", "answer_source", "prompt", "choices", "answer_clean")) %>% 
  # create index
  mutate(order_mem_index = paste(stimuli_name, pair_num, sep = "_")) %>% 
  mutate(person_trial_index = paste(run_id, order_mem_index, sep = "_")) %>% 
  # score fine-level TOM 
  mutate(Fine_TOM_result = as.numeric(response == correct_answer))

# Select only CS_across and FS_across trials
Fine_TOM_across_data <- Fine_TOM_data %>% filter(pair_type == "across") %>% 
  dplyr::select("run_id", "stimuli_type", "stimuli_name", "order_mem_index", "Fine_TOM_result")

Data exclusion: Exclude Fine_TOM trials that have reaction time < 300 ms, or > 3SD above the mean reaction time of all Fine_TOM trials

# Calculate TOM_rt upper limit: 3SD above mean 
Fine_TOM_rt_upperLim <- as.numeric(Fine_TOM_data %>% summarize(upperLim = mean(as.numeric(rt)) + 3 * sd(as.numeric(rt))))
Fine_TOM_rt_upperLim  # Print out TOM_rt upper limit: 40876.08 ms

## [1] 40876.08

# Tag all the TOM trials that are excluded by rt < 300 or > 3SD above mean
Fine_TOM_data$excluded_by_FineTOMrt <- ifelse(as.numeric(Fine_TOM_data$rt) > Fine_TOM_rt_upperLim | as.numeric(Fine_TOM_data$rt) < 300, yes = TRUE, no = FALSE)
mean(Fine_TOM_data$excluded_by_FineTOMrt == TRUE) # Print out the proportion of Fine_TOM trials being excluded based on rt criterion

## [1] 0.006377551

3.3.1 Fine-level Recency Judgment Analysis

Analysis 1

[Using all 80 within & across test-pairs; Aim to replicate results from previous experiments]

[Most important contrast: CS_across > FS_across]

Hypothesis (H1): Semantic knowledge facilitates temporal order memory. Whether the presence of event boundary impairs or facilitates temporal order memory depends on whether semantic facilitation is on coarse-level events or on fine-level events. (Interaction between narrative_type and pair_type)

• When there is semantic knowledge that helps people infer the temporal order between coarse-level events, temporal order memory of across-event pairs (CS_across pairs) will be better than when there is no semantic knowledge facilitation (FS_across and CS_within pairs).

• Similarly, when there is semantic knowledge facilitation among fine-level events, temporal order memory of within-event pairs (FS_within pairs) will be better than when there is no semantic knowledge facilitation (CS_within and FS_across pairs). - Model: Mixed Effects Logistic Regression Model - Outcome variable: Temporal Order Memory Accuracy (0/1) - Predictors: stimuli_type (CS narrative vs. FS narrative) * pair_type (across-event vs. within-event) - Random Effects: run_id (subject) + order_mem_index (different test pairs)

mean(Fine_TOM_data$Fine_TOM_result)

## [1] 0.7711735

sd(Fine_TOM_data$Fine_TOM_result)

## [1] 0.4201309

std.error(Fine_TOM_data$Fine_TOM_result)

## [1] 0.006710295

# Create Fine_TOM_data dataframe that excluded problematic Fine_TOM trials
# (Excluded 15 trials)
Fine_TOM_data_analysis <- Fine_TOM_data %>% filter(excluded_by_FineTOMrt == FALSE)
Fine_TOM_data_analysis$event_pair_type <- paste(Fine_TOM_data_analysis$stimuli_type, Fine_TOM_data_analysis$pair_type, sep = "_")

Fine_TOM_summ <- Fine_TOM_data_analysis %>% 
  group_by(run_id, event_pair_type) %>% 
  summarize(mean_FineTOM_accuracy = mean(Fine_TOM_result),
            sd_FineTOM_accuracy = sd(Fine_TOM_result))

## `summarise()` has grouped output by 'run_id'. You can override using the
## `.groups` argument.

#Fine_TOM_summ
Fine_TOM_summ_plot <- Fine_TOM_summ %>%
  ggplot(aes(x = event_pair_type, y = mean_FineTOM_accuracy, fill = event_pair_type)) +
    geom_violin() +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    geom_line(aes(group = run_id)) 

Fine_TOM_summ %>% group_by(event_pair_type) %>% summarize(total_mean = mean(mean_FineTOM_accuracy))

Fine_TOM_data_analysis$Fine_TOM_result <- as.factor(Fine_TOM_data_analysis$Fine_TOM_result)

# Contrast Coding
Fine_TOM_data_analysis_new <- Fine_TOM_data_analysis
Fine_TOM_data_analysis_new$stimuli_type <- as.factor(Fine_TOM_data_analysis_new$stimuli_type)
Fine_TOM_data_analysis_new$pair_type <- as.factor(Fine_TOM_data_analysis_new$pair_type)

contrasts(Fine_TOM_data_analysis_new$stimuli_type) <- contr.sum(2) 
contrasts(Fine_TOM_data_analysis_new$pair_type) <- contr.sum(2)
contrasts(Fine_TOM_data_analysis_new$stimuli_type)

##    [,1]
## CS    1
## FS   -1

contrasts(Fine_TOM_data_analysis_new$pair_type)

##        [,1]
## across    1
## within   -1

# Model 1: Random slope of stimuli_type * pair_type by subject, and random intercept of each_pair
TOM_accuracy_maximal_model <- glmer(Fine_TOM_result ~ stimuli_type * pair_type + (stimuli_type * pair_type |run_id) + (1|order_mem_index), data = Fine_TOM_data_analysis_new, family = binomial, control = cl1)

# Model 2: Random slopes of stimuli_type by subject, and random intercept of each_pair
TOM_accuracy_SlopeStimuliType_model <- glmer(Fine_TOM_result ~ stimuli_type * pair_type + (stimuli_type |run_id) + (1|order_mem_index), data = Fine_TOM_data_analysis_new, family = binomial, control = cl1)

# Model 3: Random slopes of pair_type by subject, and random intercept of each_pair
TOM_accuracy_SlopePairType_model <- glmer(Fine_TOM_result ~ stimuli_type * pair_type + (pair_type|run_id) + (1|order_mem_index), data = Fine_TOM_data_analysis_new, family = binomial, control = cl1)

# Model 4: Random intercepts of subject and each_pair
TOM_accuracy_Intercept_model <- glmer(Fine_TOM_result ~ stimuli_type * pair_type + (1|run_id) + (1|order_mem_index), data = Fine_TOM_data_analysis_new, family = binomial)

# Model Comparison favors Model 4 (Random intercepts of subject and each_pair)
anova(TOM_accuracy_Intercept_model, TOM_accuracy_SlopePairType_model, TOM_accuracy_SlopeStimuliType_model, TOM_accuracy_maximal_model)

TOM_accuracy_model <- TOM_accuracy_Intercept_model
summ(TOM_accuracy_model)

## MODEL INFO:
## Observations: 3895
## Dependent Variable: Fine_TOM_result
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 3735.51, BIC = 3773.12
## Pseudo-R² (fixed effects) = 0.08
## Pseudo-R² (total) = 0.28 
## 
## FIXED EFFECTS:
## -------------------------------------------------------------
##                                   Est.   S.E.   z val.      p
## ------------------------------ ------- ------ -------- ------
## (Intercept)                       1.51   0.12    12.17   0.00
## stimuli_type1                    -0.19   0.11    -1.73   0.08
## pair_type1                       -0.15   0.11    -1.39   0.16
## stimuli_type1:pair_type1          0.57   0.11     5.29   0.00
## -------------------------------------------------------------
## 
## RANDOM EFFECTS:
## -------------------------------------------
##       Group         Parameter    Std. Dev. 
## ----------------- ------------- -----------
##  order_mem_index   (Intercept)     0.86    
##      run_id        (Intercept)     0.43    
## -------------------------------------------
## 
## Grouping variables:
## -----------------------------------
##       Group        # groups   ICC  
## ----------------- ---------- ------
##  order_mem_index      80      0.18 
##      run_id           49      0.04 
## -----------------------------------

Anova(TOM_accuracy_model, type = "III")

# Check model assumptions
check_collinearity(TOM_accuracy_model) # Low Multicollinearity

#dotplot(ranef(TOM_accuracy_model,condVar=TRUE),scales = list(x =list(relation = 'free'))) # Check Random intercepts
TOM_accuracy_model_simulationOutput <- simulateResiduals(fittedModel = TOM_accuracy_model, plot = F)
plot(TOM_accuracy_model_simulationOutput)

emmean_TOM_summ <- emmeans(TOM_accuracy_model, specs = ~ stimuli_type * pair_type)
summary(emmean_TOM_summ, type = "response")

# Set up comparisons
CS_across = c(1,0,0,0)
FS_across = c(0,1,0,0)
CS_within = c(0,0,1,0)
FS_within = c(0,0,0,1)

# Planned Contrast - Interaction between stimuli_type and pair_type
contrast(emmean_TOM_summ, method = list("FS_within - CS_within" = FS_within - CS_within,
                                         "FS_within - FS_across" = FS_within - FS_across,
                                         "CS_across - FS_across" = CS_across - FS_across,
                                         "CS_across - CS_within" = CS_across - CS_within), 
         adjust = "none")

##  contrast              estimate    SE  df z.ratio p.value
##  FS_within - CS_within    1.506 0.307 Inf   4.901  <.0001
##  FS_within - FS_across    1.433 0.308 Inf   4.656  <.0001
##  CS_across - FS_across    0.764 0.299 Inf   2.558  0.0105
##  CS_across - CS_within    0.837 0.298 Inf   2.809  0.0050
## 
## Results are given on the log odds ratio (not the response) scale.

emmeans(TOM_accuracy_model, pairwise ~ stimuli_type * pair_type, adjust = "none")

## $emmeans
##  stimuli_type pair_type emmean    SE  df asymp.LCL asymp.UCL
##  CS           across     1.747 0.223 Inf     1.309      2.19
##  FS           across     0.983 0.217 Inf     0.558      1.41
##  CS           within     0.910 0.216 Inf     0.486      1.33
##  FS           within     2.416 0.236 Inf     1.954      2.88
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df z.ratio p.value
##  CS across - FS across   0.7643 0.299 Inf   2.558  0.0105
##  CS across - CS within   0.8375 0.298 Inf   2.809  0.0050
##  CS across - FS within  -0.6689 0.312 Inf  -2.147  0.0318
##  FS across - CS within   0.0732 0.294 Inf   0.249  0.8030
##  FS across - FS within  -1.4332 0.308 Inf  -4.656  <.0001
##  CS within - FS within  -1.5064 0.307 Inf  -4.901  <.0001
## 
## Results are given on the log odds ratio (not the response) scale.

# Plot contrasts - interaction between narrative_type and pair_type
plot_model(TOM_accuracy_model, type = "int", c("pair_type", "stimuli_type"), colors = c("#EE8227", "#79ADD6"), dot.size = 3.5, line.size = 1.1, axis.labels = c("CS Narrative", "FS Narrative")) + 
  scale_y_continuous(limits = c(0.5, 1), breaks = c(0.5, 0.6, 0.7, 0.8, 0.9, 1), labels = scales::percent) +
  theme_classic() +
  theme(aspect.ratio = 1.1, axis.text = element_text(size = 12)) +
  #aes(color=group) +
  #font_size(axis_title.x = 13, axis_title.y = 13,labels.x = 14, labels.y = 14) + 
  labs(x = "Narrative Type", y = "Estimated Recency Judgment Accuracy", title = NULL, size = 2) 

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Select Only "across" Trials in Fine_TOM_data
Fine_TOM_data_cleaned <- Fine_TOM_data %>% 
  filter(pair_type == "across")

## Fine_TOM_Conf data
Fine_TOM_Conf_data <- main_data %>% filter(task == "TOM_Conf")
Fine_TOM_Conf_data <- Fine_TOM_Conf_data %>% 
  # select only useful columns
  dplyr::select(-c("SourceMem_question_index", "pair_within_index", "coarse_count", "SourceMem_correct_answer", "Source_Sentence_index", "SourceMem_correct", "answer_source", "prompt", "choices", "answer_clean")) %>% 
  # create index
  mutate(order_mem_index = paste(stimuli_name, pair_num, sep = "_")) %>% 
  mutate(person_trial_index = paste(run_id, order_mem_index, sep = "_")) 

Fine_TOM_Conf_data_cleaned <- Fine_TOM_Conf_data %>% 
  # only select across trials! 
  filter(pair_type == "across") %>% 
  dplyr::select("person_trial_index", "response")  
colnames(Fine_TOM_Conf_data_cleaned)[colnames(Fine_TOM_Conf_data_cleaned) == "response"] <- "Fine_TOM_Conf"

# MERGE: Fine_TOM (total) + Fine_TOM_Conf
Fine_TOM_complete_data_cleaned <- merge(Fine_TOM_data_cleaned, Fine_TOM_Conf_data_cleaned, by = "person_trial_index") 

3.4 Merge Data

Merge Fine_TOM, Fine_TOM_Conf, Coarse_TOM, Coarse_TOM_Conf, Source_Mem trials according to the common “person_trial_index”.

# Step 1: MERGE Fine_TOM_data_cleaned with source_mem_data_cleaned
Fine_Source_merged <- merge(Fine_TOM_complete_data_cleaned, source_mem_data_cleaned, by = "person_trial_index") 

# Step2: MERGE Fine_Source_merged with Coarse_TOM_data_cleaned
Fine_Source_Coarse_merged <- merge(Fine_Source_merged, Coarse_TOM_complete_data_cleaned, by = "person_trial_index") 

FSC_data <- Fine_Source_Coarse_merged %>% dplyr::select ("run_id", "stimuli_order", "stimuli_type", "pair_type", "stimuli_name",  "pair_first", "pair_second", "Coarse_pair_first", "Coarse_pair_second", "order_mem_index", "person_trial_index", "Fine_TOM_result", "Fine_TOM_Conf", "source_mem_result_s1", "source_mem_result_s2", "answer_source_s1", "answer_source_s2", "correct_source_s1", "correct_source_s2", "Coarse_TOM_result", "Coarse_TOM_Conf", "excluded_by_FineTOMrt", "excluded_by_Coarse_TOM_rt", "excluded_by_sourceMem1_rt", "excluded_by_sourceMem2_rt") 

4. Analysis of Across Trials only (40 pairs)

# Exclude Problematic Trials according to rt (19 out of 560, with 541 left)
FSC_data <- FSC_data %>% filter(excluded_by_Coarse_TOM_rt == FALSE & excluded_by_FineTOMrt == FALSE & excluded_by_sourceMem1_rt == FALSE & excluded_by_sourceMem2_rt == FALSE)

# Graph the distribution of Fine_TOM Accuracy across event_pair_types
Fine_TOM_data_dist <- ggplot(FSC_data, aes(x = as.factor(Fine_TOM_result), fill = stimuli_type)) + 
  geom_bar() +
  labs(title = "Distribution of Fine_TOM Accuracy grouped by event pair types", x = "Fine_TOM Accuracy", y = "Count")
#Fine_TOM_data_dist

# Graph the distribution of Coarse_TOM Accuracy across event_pair_types
Coarse_TOM_data_dist <- ggplot(FSC_data, aes(x = as.factor(Coarse_TOM_result), fill = stimuli_type)) + 
  geom_bar() +
  labs(title = "Distribution of Coarse_TOM Accuracy grouped by event pair types", x = "Fine_TOM Accuracy", y = "Count")
#Coarse_TOM_data_dist

FSC_data %>% 
  group_by(stimuli_type) %>% 
  summarize(mean_coarse_accuracy = mean(Coarse_TOM_result),
            count_coarse_0 = sum(Coarse_TOM_result == 0),
            count_coarse_1 = sum(Coarse_TOM_result == 1))

# Graph the distribution of Temporal Order Memory Trial Accuracy across event_pair_types
SourceMem_s1_data_dist <- ggplot(FSC_data, aes(x = as.factor(source_mem_result_s1), fill = stimuli_type)) + 
  geom_bar() +
  labs(title = "Distribution of SourceMem_s1 Accuracy grouped by event pair types", x = "SourceMem_s1 Accuracy", y = "Count")

# Graph the distribution of Temporal Order Memory Trial Accuracy across event_pair_types
SourceMem_s2_data_dist <- ggplot(FSC_data, aes(x = as.factor(source_mem_result_s2), fill = stimuli_type)) + 
  geom_bar() +
  labs(title = "Distribution of SourceMem_s2 Accuracy grouped by event pair types", x = "SourceMem_s2 Accuracy", y = "Count")
#SourceMem_s2_data_dist

FSC_data %>% 
  group_by(stimuli_type) %>% 
  summarize(mean_s1_accuracy = mean(source_mem_result_s1),
            mean_s2_accuracy = mean(source_mem_result_s2),
            count_s1_0 = sum(source_mem_result_s1 == 0),
            count_s2_0 = sum(source_mem_result_s2 == 0),
            count_s1_1 = sum(source_mem_result_s1 == 1),
            count_s2_1 = sum(source_mem_result_s2 == 1))

4.1 Analysis of CS/FS difference on Fine_TOM results

Analysis 2

Note: We did not separately report this analysis in the manuscript, because this reported the same results as the CS_across > FS_across contrast in the fine-level recency judgment analysis above.

[Using 40 across-event fine-level test-pairs; 20 CS_across, 20 FS_across]

Hypothesis (H2): Coarse-level semantic knowledge facilitates fine-level temporal order memory. When there is semantic knowledge that helps people infer the temporal order between coarse-level events, temporal order memory of across-event pairs (CS_across pairs) will be better than when there is no semantic knowledge facilitation (FS_across pairs).

• Model: Mixed Effects Logistic Regression Model
• Outcome variable: Fine-level Temporal Order Memory Accuracy (0/1)
• Predictors: stimuli_type (CS narrative vs. FS narrative)
• Grouping variables of Random Effects: run_id (subject) + order_mem_index (different test pairs)

Fine_TOM_summ <- FSC_data %>% 
  group_by(run_id, stimuli_type) %>% 
  summarize(mean_FineTOM_accuracy = mean(Fine_TOM_result),
            sd_FineTOM_accuracy = sd(Fine_TOM_result))

## `summarise()` has grouped output by 'run_id'. You can override using the
## `.groups` argument.

Fine_TOM_summ

Fine_TOM_summ %>%
  ggplot(aes(x = stimuli_type, y = mean_FineTOM_accuracy, fill = stimuli_type)) +
    geom_violin() +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    geom_line(aes(group = run_id)) 

Fine_TOM_summ %>% group_by(stimuli_type) %>% summarize(total_mean = mean(mean_FineTOM_accuracy))

# Make binary variables
FSC_data$Fine_TOM_result_binary <- as.factor(FSC_data$Fine_TOM_result)
FSC_data$Coarse_TOM_result_binary <- as.factor(FSC_data$Coarse_TOM_result)
FSC_data$stimuli_type <- as.factor(FSC_data$stimuli_type)

# Contrast Coding
contrasts(FSC_data$stimuli_type) <- contr.sum(2)
contrasts(FSC_data$stimuli_type)

##    [,1]
## CS    1
## FS   -1

# Model 1: With stimuli_type random slope on both subject and event pair
Fine_TOM_model_maximal <- glmer(Fine_TOM_result_binary ~ 1 + stimuli_type + (1 + stimuli_type | run_id) + (1 + stimuli_type | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 2: With stimuli_type random slope on only subject
Fine_TOM_model_subjSlope <- glmer(Fine_TOM_result_binary ~ 1 + stimuli_type + (1 + stimuli_type | run_id) + (1 | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 3: With stimuli_type random slope on only event pair
Fine_TOM_model_pairSlope <- glmer(Fine_TOM_result_binary ~ 1 + stimuli_type + (1 | run_id) + (1 + stimuli_type | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 4: With only random intercepts 
Fine_TOM_model_intercept <- glmer(Fine_TOM_result_binary ~ 1 + stimuli_type + (1 | run_id) + (1 | order_mem_index), data = FSC_data, family = binomial)

# Model selection favors Model 4
anova(Fine_TOM_model_intercept, Fine_TOM_model_subjSlope, Fine_TOM_model_pairSlope, Fine_TOM_model_maximal)

# Fine_TOM_model Output:
Fine_TOM_model <- Fine_TOM_model_intercept
summ(Fine_TOM_model)

## MODEL INFO:
## Observations: 1915
## Dependent Variable: Fine_TOM_result_binary
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 1930.01, BIC = 1952.24
## Pseudo-R² (fixed effects) = 0.04
## Pseudo-R² (total) = 0.25 
## 
## FIXED EFFECTS:
## -------------------------------------------------
##                       Est.   S.E.   z val.      p
## ------------------- ------ ------ -------- ------
## (Intercept)           1.38   0.17     8.36   0.00
## stimuli_type1         0.39   0.15     2.58   0.01
## -------------------------------------------------
## 
## RANDOM EFFECTS:
## -------------------------------------------
##       Group         Parameter    Std. Dev. 
## ----------------- ------------- -----------
##      run_id        (Intercept)     0.42    
##  order_mem_index   (Intercept)     0.88    
## -------------------------------------------
## 
## Grouping variables:
## -----------------------------------
##       Group        # groups   ICC  
## ----------------- ---------- ------
##      run_id           49      0.04 
##  order_mem_index      40      0.18 
## -----------------------------------

Anova(Fine_TOM_model, type = "III")

4.2 Analysis of CS/FS difference on Coarse_TOM results

Analysis 3

[Using 40 coarse-level test-pairs; 20 CS_across, 20 FS_across]

Hypothesis (H3): Coarse-level semantic knowledge facilitates coarse-level temporal order memory. When there is semantic knowledge that helps people infer the temporal order between coarse-level events, temporal order memory of across-event pairs (CS pairs) will be better than when there is no semantic knowledge facilitation (FS pairs).

• Model: Mixed Effects Logistic Regression Model
• Outcome variable: Fine-level Temporal Order Memory Accuracy (0/1)
• Predictors: stimuli_type (CS narrative vs. FS narrative)
• Grouping variables of Random Effects: run_id (subject) + order_mem_index (different test pairs)

mean(FSC_data$Coarse_TOM_result)

## [1] 0.8496084

sd(FSC_data$Coarse_TOM_result)

## [1] 0.3575483

std.error(FSC_data$Coarse_TOM_result)

## [1] 0.008170531

Coarse_TOM_summ <- FSC_data %>% 
  group_by(run_id, stimuli_type) %>% 
  summarize(mean_CoarseTOM_accuracy = mean(Coarse_TOM_result),
            sd_CoarseTOM_accuracy = sd(Coarse_TOM_result))

## `summarise()` has grouped output by 'run_id'. You can override using the
## `.groups` argument.

# Coarse_TOM_summ
Coarse_TOM_summ %>%
  ggplot(aes(x = stimuli_type, y = mean_CoarseTOM_accuracy, fill = stimuli_type)) +
    geom_violin() +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    geom_line(aes(group = run_id)) 

Coarse_TOM_summ %>% group_by(stimuli_type) %>% summarize(total_mean = mean(mean_CoarseTOM_accuracy))

# Model 1: With stimuli_type random slope on both subject and event pair
Coarse_TOM_model_maximal <- glmer(Coarse_TOM_result_binary ~ 1 + stimuli_type + (1 + stimuli_type | run_id) + (1 + stimuli_type | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 2: With stimuli_type random slope on only subject
Coarse_TOM_model_subjSlope <- glmer(Coarse_TOM_result_binary ~ 1 + stimuli_type + (1 + stimuli_type | run_id) + (1 | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 3: With stimuli_type random slope on only event pair
Coarse_TOM_model_pairSlope <- glmer(Coarse_TOM_result_binary ~ 1 + stimuli_type + (1 | run_id) + (1 + stimuli_type| order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model 4: With only random intercepts 
Coarse_TOM_model_intercept <- glmer(Coarse_TOM_result_binary ~ 1 + stimuli_type + (1 | run_id) + (1 | order_mem_index), data = FSC_data, family = binomial, control = cl1)

# Model selection favors Model 2
anova(Coarse_TOM_model_intercept, Coarse_TOM_model_subjSlope, Coarse_TOM_model_pairSlope, Coarse_TOM_model_maximal)

anova(Coarse_TOM_model_subjSlope, Coarse_TOM_model_maximal)

# Coarse_TOM_model Output:
Coarse_TOM_model <- Coarse_TOM_model_subjSlope
summ(Coarse_TOM_model)

## MODEL INFO:
## Observations: 1915
## Dependent Variable: Coarse_TOM_result_binary
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 1401.18, BIC = 1434.53
## Pseudo-R² (fixed effects) = 0.23
## Pseudo-R² (total) = 0.47 
## 
## FIXED EFFECTS:
## -------------------------------------------------
##                       Est.   S.E.   z val.      p
## ------------------- ------ ------ -------- ------
## (Intercept)           2.52   0.20    12.36   0.00
## stimuli_type1         1.21   0.18     6.85   0.00
## -------------------------------------------------
## 
## RANDOM EFFECTS:
## ---------------------------------------------
##       Group          Parameter     Std. Dev. 
## ----------------- --------------- -----------
##      run_id         (Intercept)      0.80    
##      run_id        stimuli_type1     0.40    
##  order_mem_index    (Intercept)      0.84    
## ---------------------------------------------
## 
## Grouping variables:
## -----------------------------------
##       Group        # groups   ICC  
## ----------------- ---------- ------
##      run_id           49      0.14 
##  order_mem_index      40      0.15 
## -----------------------------------

Anova(Coarse_TOM_model, type = "III")

# Check model assumptions
Coarse_TOM_model_simulationOutput <- simulateResiduals(fittedModel = Coarse_TOM_model, plot = F)
plot(Coarse_TOM_model_simulationOutput)

plot_model(Coarse_TOM_model, type = "emm", terms = c("stimuli_type"), colors = c("red"), dot.size = 3.5, line.size = 1.1) +
  scale_y_continuous(limits = c(0.5, 1), breaks = c(0.5, 0.6, 0.7, 0.8, 0.9, 1), labels = scales::percent) +
  theme_classic() + 
  theme(aspect.ratio = 1.1, axis.text = element_text(size = 12)) +
  labs(x = "Narrative Type", y = "Estimated Coarse-level Recency Judgement Accuracy", title = NULL, size = 2)

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

5. Analysis of FS_across trials only (20 pairs)

5.1 Fine-level FS_across Accuracy

Analysis 4

[Using 20 FS_across fine-level test pairs, 20 corresponding coarse-level test-pairs, and 40 source memory trials of sentences from the 20 fine-level test pairs]

Question: Is FS_across Fine_TOM accuracy predicted by Coarse_TOM accuracy, controlling for source_both_correct?

Hypothesis: The accuracy of a given fine-level temporal order memory trial will be predicted by whether its corresponding coarse-level order memory trial is correct, controlling for the source memory of the two fine-level event sentences.

• Model: Mixed Effects Logistic Regression Model
• Outcome variable: Fine-level Temporal Order Memory Accuracy (0/1)
• Predictors: Coarse-level Temporal Order Memory Accuracy (0/1) and whether the source memory of the two fine-level event sentences are both correct (0/1)
• Grouping variables of Random Effects: run_id (subject) + order_mem_index (different test pairs)

# Create source_both_correct variable; 1 if both s1 and s2 is correct, 0 otherwise
FSC_data <- FSC_data %>% 
  mutate(source_both_correct = 
           ifelse(source_mem_result_s1 == 1 & source_mem_result_s2 == 1, 1, 0))

FSC_data$source_both_correct <- as.factor(FSC_data$source_both_correct)

#table(FSC_data$stimuli_type, FSC_data$source_mem_result_s1)
#table(FSC_data$stimuli_type, FSC_data$source_mem_result_s2)

table(FSC_data$stimuli_type, FSC_data$source_both_correct)

##     
##        0   1
##   CS 218 736
##   FS 207 754

table(FSC_data$stimuli_type, FSC_data$Coarse_TOM_result)

##     
##        0   1
##   CS  55 899
##   FS 233 728

FS_FSC_data <- FSC_data %>% filter(stimuli_type == "FS")
CS_FSC_data <- FSC_data %>% filter(stimuli_type == "CS")

CS_FSC_data$SC_group <- paste("S", CS_FSC_data$source_both_correct, "C", CS_FSC_data$Coarse_TOM_result, sep = "_")
FS_FSC_data$SC_group <- paste("S", FS_FSC_data$source_both_correct, "C", FS_FSC_data$Coarse_TOM_result, sep = "_")

FS_FSC_summ <- FS_FSC_data %>% 
  group_by(run_id, SC_group) %>% 
  summarize(mean_FineTOM_accuracy = mean(Fine_TOM_result),
            sd_FineTOM_accuracy = sd(Fine_TOM_result),
            n = length(Fine_TOM_result))

## `summarise()` has grouped output by 'run_id'. You can override using the
## `.groups` argument.

#FS_FSC_summ
FS_FSC_group_summ_plot <- FS_FSC_summ %>%
  ggplot(aes(x = SC_group, y = mean_FineTOM_accuracy, fill = SC_group)) +
    geom_violin() +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    geom_line(aes(group = run_id)) 

FS_FSC_summ %>% group_by(SC_group) %>% summarize(total_mean = mean(mean_FineTOM_accuracy))

FS_FSC_group_summ <- FS_FSC_data %>% 
  group_by(SC_group, person_trial_index) %>% 
  summarize(mean_FineTOM_accuracy = mean(Fine_TOM_result),
            sd_FineTOM_accuracy = sd(Fine_TOM_result),
            n = length(Fine_TOM_result)) 

## `summarise()` has grouped output by 'SC_group'. You can override using the
## `.groups` argument.

table(FS_FSC_data$SC_group, FS_FSC_data$Fine_TOM_result)

##          
##             0   1
##   S_0_C_0  32  31
##   S_0_C_1  72  72
##   S_1_C_0 113  57
##   S_1_C_1  77 507

# Effect Coding
FS_FSC_data_new <- FS_FSC_data
FS_FSC_data_new$Coarse_TOM_result_binary <- factor(FS_FSC_data_new$Coarse_TOM_result_binary, levels = c(1, 0))
FS_FSC_data_new$source_both_correct <- factor(FS_FSC_data_new$source_both_correct, levels = c(1, 0))

contrasts(FS_FSC_data_new$Coarse_TOM_result_binary) <- contr.sum(2) 
contrasts(FS_FSC_data_new$source_both_correct) <- contr.sum(2)
contrasts(FS_FSC_data_new$Coarse_TOM_result_binary)

##   [,1]
## 1    1
## 0   -1

contrasts(FS_FSC_data_new$source_both_correct)

##   [,1]
## 1    1
## 0   -1

# FS across: if source & coarse jointly predict fine

# Model 1: Maximal Model with interaction
Fine_TOM_FSa_model_int_maximal <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (Coarse_TOM_result_binary * source_both_correct| run_id) + (Coarse_TOM_result_binary * source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model 2: Maximal without interaction
Fine_TOM_FSa_model_noint_maximal <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + (Coarse_TOM_result_binary + source_both_correct| run_id) + (Coarse_TOM_result_binary + source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model Comparison favors Model with Interaction
anova(Fine_TOM_FSa_model_noint_maximal, Fine_TOM_FSa_model_int_maximal)

# Model a: 
Fine_TOM_FSa_model_int_a <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + Coarse_TOM_result_binary + source_both_correct| run_id) + (1 + Coarse_TOM_result_binary + source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

anova(Fine_TOM_FSa_model_int_a, Fine_TOM_FSa_model_int_maximal)

# Model b: 
Fine_TOM_FSa_model_int_b <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + Coarse_TOM_result_binary | run_id) + (1 + Coarse_TOM_result_binary + source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model c: 
Fine_TOM_FSa_model_int_c <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + source_both_correct| run_id) + (1 + Coarse_TOM_result_binary + source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model d: 
Fine_TOM_FSa_model_int_d <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 | run_id) + (1 + Coarse_TOM_result_binary + source_both_correct|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

anova(Fine_TOM_FSa_model_int_d, Fine_TOM_FSa_model_int_c, Fine_TOM_FSa_model_int_b, Fine_TOM_FSa_model_int_a)

# Model e: 
Fine_TOM_FSa_model_int_e <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + Coarse_TOM_result_binary | run_id) + (1 + Coarse_TOM_result_binary|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model f: 
Fine_TOM_FSa_model_int_f <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + Coarse_TOM_result_binary | run_id) + (1 + source_both_correct | order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model g: 
Fine_TOM_FSa_model_int_g <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 + Coarse_TOM_result_binary | run_id) + (1 | order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model h: intercept only
Fine_TOM_FSa_model_int_h <- glmer(Fine_TOM_result_binary ~ 1 + Coarse_TOM_result_binary + source_both_correct + Coarse_TOM_result_binary:source_both_correct + (1 | run_id) + (1 |order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Select Model g
anova(Fine_TOM_FSa_model_int_h, Fine_TOM_FSa_model_int_g, Fine_TOM_FSa_model_int_f, Fine_TOM_FSa_model_int_e, Fine_TOM_FSa_model_int_b)

Fine_TOM_FSa_model <- Fine_TOM_FSa_model_int_g
summ(Fine_TOM_FSa_model)

## MODEL INFO:
## Observations: 961
## Dependent Variable: Fine_TOM_result_binary
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 917.82, BIC = 956.77
## Pseudo-R² (fixed effects) = 0.24
## Pseudo-R² (total) = 0.39 
## 
## FIXED EFFECTS:
## ----------------------------------------------------------------------------------
##                                                        Est.   S.E.   z val.      p
## ---------------------------------------------------- ------ ------ -------- ------
## (Intercept)                                            0.41   0.21     2.02   0.04
## Coarse_TOM_result_binary1                              0.72   0.12     5.85   0.00
## source_both_correct1                                   0.20   0.13     1.60   0.11
## Coarse_TOM_result_binary1:source_both_correct1         0.67   0.10     6.52   0.00
## ----------------------------------------------------------------------------------
## 
## RANDOM EFFECTS:
## ---------------------------------------------------------
##       Group                Parameter           Std. Dev. 
## ----------------- --------------------------- -----------
##      run_id               (Intercept)            0.00    
##      run_id        Coarse_TOM_result_binary1     0.49    
##  order_mem_index          (Intercept)            0.77    
## ---------------------------------------------------------
## 
## Grouping variables:
## -----------------------------------
##       Group        # groups   ICC  
## ----------------- ---------- ------
##      run_id           49      0.00 
##  order_mem_index      20      0.15 
## -----------------------------------

Anova(Fine_TOM_FSa_model, type = "III")

# Check model assumptions
check_collinearity(Fine_TOM_FSa_model) # Low Multicollinearity

Fine_TOM_FSa_model_simulationOutput <- simulateResiduals(fittedModel = Fine_TOM_FSa_model, plot = F)
plot(Fine_TOM_FSa_model_simulationOutput)

# Probe the interaction with pairwise comparisons - Simple effects
emmeans(Fine_TOM_FSa_model, pairwise ~ Coarse_TOM_result_binary|source_both_correct, adjust = "bonferroni")

## $emmeans
## source_both_correct = 1:
##  Coarse_TOM_result_binary emmean    SE  df asymp.LCL asymp.UCL
##  1                         2.008 0.231 Inf     1.556     2.460
##  0                        -0.778 0.265 Inf    -1.297    -0.259
## 
## source_both_correct = 0:
##  Coarse_TOM_result_binary emmean    SE  df asymp.LCL asymp.UCL
##  1                         0.268 0.291 Inf    -0.302     0.839
##  0                         0.157 0.356 Inf    -0.541     0.855
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
## source_both_correct = 1:
##  contrast                                              estimate    SE  df
##  Coarse_TOM_result_binary1 - Coarse_TOM_result_binary0    2.787 0.268 Inf
##  z.ratio p.value
##   10.392  <.0001
## 
## source_both_correct = 0:
##  contrast                                              estimate    SE  df
##  Coarse_TOM_result_binary1 - Coarse_TOM_result_binary0    0.111 0.368 Inf
##  z.ratio p.value
##    0.303  0.7621
## 
## Results are given on the log odds ratio (not the response) scale.

# Plot contrasts - interaction between narrative_type and pair_type
plot_model(Fine_TOM_FSa_model, type = "int", colors = c("#AD5691", "#EE3224"), dot.size = 3.5, line.size = 1.1, legend.title = "Source Memory (Both Correct)", legend.labels = c("False", "True")) + 
  theme_classic() + 
  theme(aspect.ratio = 1.1, axis.text = element_text(size = 12)) +
  coord_cartesian(ylim = c(0.2, 1)) +
  #aes(color=group) +
  #font_size(axis_title.x = 13, axis_title.y = 13,labels.x = 14, labels.y = 14) + 
  labs(x = "Coarse-level Recency Judgment Accuracy", y = "Estimated Recency Judgment Accuracy", title = NULL, size = 2) 

5.2 Fine-level FS_across Confidence

Analysis 5

[Using 20 FS_across fine-level test pairs, 20 corresponding coarse-level test-pairs, and 40 source memory trials of sentences from the 20 fine-level test pairs]

Question: Is FS_across Fine_TOM Confidence predicted by Coarse_TOM Confidence, controlling for Fine_TOM_result, Coarse_TOM_result, and source_both_correct?

Hypothesis: We hypothesize that the confidence of a given fine-level recency judgment trial (binary, high/low) will be more strongly predicted by the confidence of its corresponding coarse-level recency judgment trial (binary, high/low) when the fine-level recency judgment is correct, controlling for the accuracy of coarse-level recency judgment (0/1) and the source memory of the two fine-level event sentences (0/1).

• Model: Mixed Effects Logistic Regression Model
• Outcome variable: Fine-level Temporal Order Memory Confidence (High/Low Confidence)
• Predictors: Coarse-level Temporal Order Memory Confidence (High/Low Confidence), Fine-level Temporal Order Memory Accuracy (0/1), Coarse-level Temporal Order Memory Accuracy (0/1), and whether the source memory of the two fine-level event sentences are both correct (0/1)
• Grouping variables of Random Effects: run_id (subject) + order_mem_index (different test pairs)

# Graph the distribution of Fine-level Temporal Order Memory Confidence across event_pair_types
Fine_TOM_Conf_data_dist <- ggplot(FS_FSC_data, aes(x = as.numeric(Fine_TOM_Conf), fill = stimuli_type)) + 
  geom_histogram() +   
  labs(title = "Distribution of Fine_TOM Confidence grouped by event pair types", x = "Fine_TOM Confidence", y = "Count")
#Fine_TOM_Conf_data_dist

FS_FSC_data$Fine_TOM_Conf <- as.numeric(FS_FSC_data$Fine_TOM_Conf)
FS_FSC_data$Coarse_TOM_Conf <- as.numeric(FS_FSC_data$Coarse_TOM_Conf)

mean(FS_FSC_data$Fine_TOM_Conf)

## [1] 69.26431

sd(FS_FSC_data$Fine_TOM_Conf)

## [1] 26.95928

std.error(FS_FSC_data$Fine_TOM_Conf)

## [1] 0.8696543

# Graph the distribution of Coarse-level Temporal Order Memory Confidence across event_pair_types
Coarse_TOM_Conf_data_dist <- ggplot(FS_FSC_data, aes(x = as.numeric(Coarse_TOM_Conf), fill = stimuli_type)) + 
  geom_histogram() +   
  labs(title = "Distribution of Coarse_TOM Confidence grouped by event pair types", x = "Coarse_TOM Confidence", y = "Count")
#Coarse_TOM_Conf_data_dist

mean(FS_FSC_data$Coarse_TOM_Conf)

## [1] 72.55359

sd(FS_FSC_data$Coarse_TOM_Conf)

## [1] 26.76338

std.error(FS_FSC_data$Coarse_TOM_Conf)

## [1] 0.8633349

# Since 100 is the mode, we use 80 as the cutoff for converting TOM_Conf into a binary variable 
# Binary TOM Confidence: "High Confidence" (Confidence > 80) vs. "Low Confidence" (Confidence < 80)
FS_FSC_data$Fine_TOM_ConfGroup <- ifelse(test = FS_FSC_data$Fine_TOM_Conf > 80, yes = "High Confidence", no = "Low Confidence")
FS_FSC_data$Coarse_TOM_ConfGroup <- ifelse(test = FS_FSC_data$Coarse_TOM_Conf > 80, yes = "High Confidence", no = "Low Confidence")

# If use 80 as cutoff, 42% FS Fine_level TOM trials are in the High Confidence group, 47% FS Coarse-level TOM trials are in the High Confidence group
mean(FS_FSC_data$Fine_TOM_ConfGroup == "High Confidence") 

## [1] 0.4203954

mean(FS_FSC_data$Coarse_TOM_ConfGroup == "High Confidence") 

## [1] 0.4651405

table(FS_FSC_data$Fine_TOM_ConfGroup)

## 
## High Confidence  Low Confidence 
##             404             557

table(FS_FSC_data$Coarse_TOM_ConfGroup)

## 
## High Confidence  Low Confidence 
##             447             514

FS_FSC_data <- FS_FSC_data %>% mutate(Fine_TOM_ConfGroup_dumm = dplyr::recode(Fine_TOM_ConfGroup,
                                             "Low Confidence" = 0,
                                             "High Confidence" = 1),
                                      Coarse_TOM_ConfGroup_dumm = dplyr::recode(Coarse_TOM_ConfGroup,
                                             "Low Confidence" = 0,
                                             "High Confidence" = 1))

# Effect Coding
FS_FSC_data_new <- FS_FSC_data
FS_FSC_data_new$Fine_TOM_result_binary <- factor(FS_FSC_data_new$Fine_TOM_result_binary, levels = c(1, 0))
FS_FSC_data_new$Coarse_TOM_result_binary <- factor(FS_FSC_data_new$Coarse_TOM_result_binary, levels = c(1, 0))
FS_FSC_data_new$Coarse_TOM_ConfGroup_dumm <- factor(FS_FSC_data_new$Coarse_TOM_ConfGroup_dumm, levels = c(1, 0))
FS_FSC_data_new$source_both_correct <- factor(FS_FSC_data_new$source_both_correct, levels = c(1, 0))

contrasts(FS_FSC_data_new$Fine_TOM_result_binary) <- contr.sum(2)
contrasts(FS_FSC_data_new$Coarse_TOM_result_binary) <- contr.sum(2)
contrasts(FS_FSC_data_new$Coarse_TOM_ConfGroup_dumm) <- contr.sum(2) 
contrasts(FS_FSC_data_new$source_both_correct) <- contr.sum(2)

contrasts(FS_FSC_data_new$Fine_TOM_result_binary)

##   [,1]
## 1    1
## 0   -1

contrasts(FS_FSC_data_new$Coarse_TOM_result_binary) 

##   [,1]
## 1    1
## 0   -1

contrasts(FS_FSC_data_new$Coarse_TOM_ConfGroup_dumm) 

##   [,1]
## 1    1
## 0   -1

contrasts(FS_FSC_data_new$source_both_correct) 

##   [,1]
## 1    1
## 0   -1

FS_FSC_data_new$Fine_TOM_ConfGroup_dumm <- as.factor(FS_FSC_data_new$Fine_TOM_ConfGroup_dumm)
contrasts(FS_FSC_data_new$Fine_TOM_ConfGroup_dumm)

##   1
## 0 0
## 1 1

# Model a: 1 Model with 4-way interaction
Fine_TOM_Conf_FSa_model_a <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * source_both_correct * Coarse_TOM_ConfGroup_dumm * Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

# Model b1-b4: 4 Models with 3-way interactions
Fine_TOM_Conf_FSa_model_b1 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * source_both_correct * Coarse_TOM_ConfGroup_dumm + Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_b2 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * source_both_correct * Coarse_TOM_result_binary + Coarse_TOM_ConfGroup_dumm + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_b3 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * Coarse_TOM_ConfGroup_dumm * Coarse_TOM_result_binary + source_both_correct + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_b4 <- glmer(Fine_TOM_ConfGroup_dumm ~ source_both_correct * Coarse_TOM_ConfGroup_dumm * Coarse_TOM_result_binary + Fine_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model c1-c6: 6 Models with 2-way interactions
Fine_TOM_Conf_FSa_model_c1 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * source_both_correct + Coarse_TOM_ConfGroup_dumm + Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_c2 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * Coarse_TOM_ConfGroup_dumm + source_both_correct + Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_c3 <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary * Coarse_TOM_result_binary + Coarse_TOM_ConfGroup_dumm + source_both_correct + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_c4 <- glmer(Fine_TOM_ConfGroup_dumm ~ source_both_correct * Coarse_TOM_ConfGroup_dumm + Fine_TOM_result_binary + Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_c5 <- glmer(Fine_TOM_ConfGroup_dumm ~ source_both_correct * Coarse_TOM_result_binary + Fine_TOM_result_binary + Coarse_TOM_ConfGroup_dumm + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

Fine_TOM_Conf_FSa_model_c6 <- glmer(Fine_TOM_ConfGroup_dumm ~ Coarse_TOM_ConfGroup_dumm * Coarse_TOM_result_binary + Fine_TOM_result_binary + source_both_correct + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model d: 1 Model with no interactions
Fine_TOM_Conf_FSa_model_d <- glmer(Fine_TOM_ConfGroup_dumm ~ Fine_TOM_result_binary + source_both_correct + Coarse_TOM_ConfGroup_dumm + Coarse_TOM_result_binary + (1 | run_id) + (1|order_mem_index), data = FS_FSC_data_new, family = binomial, control = cl1)

# Model comparison: 
anova(Fine_TOM_Conf_FSa_model_d, Fine_TOM_Conf_FSa_model_c1, Fine_TOM_Conf_FSa_model_c2, Fine_TOM_Conf_FSa_model_c3, Fine_TOM_Conf_FSa_model_c4, Fine_TOM_Conf_FSa_model_c5, Fine_TOM_Conf_FSa_model_c6, Fine_TOM_Conf_FSa_model_b1, Fine_TOM_Conf_FSa_model_b2, Fine_TOM_Conf_FSa_model_b3, Fine_TOM_Conf_FSa_model_b4, Fine_TOM_Conf_FSa_model_a)

anova(Fine_TOM_Conf_FSa_model_d, Fine_TOM_Conf_FSa_model_c3, Fine_TOM_Conf_FSa_model_b3, Fine_TOM_Conf_FSa_model_a)

# Model selection favors Fine_TOM_Conf_FSa_model_d
Fine_TOM_Conf_FSa_model <- Fine_TOM_Conf_FSa_model_d
summ(Fine_TOM_Conf_FSa_model, digits = getOption("jtools-digits", 4))

## MODEL INFO:
## Observations: 961
## Dependent Variable: Fine_TOM_ConfGroup_dumm
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit 
## 
## MODEL FIT:
## AIC = 996.2038, BIC = 1030.2796
## Pseudo-R² (fixed effects) = 0.2052
## Pseudo-R² (total) = 0.4270 
## 
## FIXED EFFECTS:
## ----------------------------------------------------------------------
##                                       Est.     S.E.    z val.        p
## -------------------------------- --------- -------- --------- --------
## (Intercept)                        -0.6095   0.2306   -2.6429   0.0082
## Fine_TOM_result_binary1             0.4187   0.1106    3.7867   0.0002
## source_both_correct1               -0.0254   0.1331   -0.1905   0.8489
## Coarse_TOM_ConfGroup_dumm1          0.9110   0.0959    9.4951   0.0000
## Coarse_TOM_result_binary1           0.0415   0.1150    0.3612   0.7179
## ----------------------------------------------------------------------
## 
## RANDOM EFFECTS:
## -------------------------------------------
##       Group         Parameter    Std. Dev. 
## ----------------- ------------- -----------
##      run_id        (Intercept)    0.9617   
##  order_mem_index   (Intercept)    0.5904   
## -------------------------------------------
## 
## Grouping variables:
## -------------------------------------
##       Group        # groups    ICC   
## ----------------- ---------- --------
##      run_id           49      0.2027 
##  order_mem_index      20      0.0764 
## -------------------------------------

Anova(Fine_TOM_Conf_FSa_model, type = "III")

# Check model assumptions
check_collinearity(Fine_TOM_Conf_FSa_model) # Low Multicollinearity

Fine_TOM_Conf_FSa_model_simulationOutput <- simulateResiduals(fittedModel = Fine_TOM_Conf_FSa_model, plot = F)
plot(Fine_TOM_Conf_FSa_model_simulationOutput)

plot_model(Fine_TOM_Conf_FSa_model, type = "emm", terms = c("Coarse_TOM_ConfGroup_dumm"), colors = c("red"), dot.size = 3.5, line.size = 1.1) +
  theme_classic() + 
  theme(aspect.ratio = 1.1, axis.text = element_text(size = 12)) +
  scale_y_continuous(limits = seq(0, 1), breaks = seq(0, 1, by = 0.2), labels = scales::percent) +
  labs(x = "Coarse-level Recency Confidence (0 = low, 1 = high)", y = "Estimated Probabilities of High Confidence (>80)", title = NULL, size = 2)

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.